The effect of \(U_{[a]}^{\text{dyad}}\) on \(\phi\) are encoded using a vector of parameters \(\beta_{\phi}\). The values \(\beta_{\phi}\) are chosen such that dyads that spend less time in state \(1\) (no body contact) spend longer in body contact and grooming (\(S \in \{2, 3, 4\}\)), and are more likely to transition from “affiliative” states (\(S \in \{2, 3, 4\}\)) than to “non-affiliative” ones (\(S = 1\)).
The effect of \(U_{[a, b]}^{\text{dyad}}\) (with \(U_{[a, b]}^{\text{dyad}} = U_{[b, a]}^{\text{dyad}}\)) on \(\tau\) is similarly encoded with a vector \(\beta_{\tau}\).
Given the population-level parameters \(\alpha\), \(\beta_{\phi}\), and \(\beta_{\tau}\), we generate the average holding times \(\theta_{[a, b, k]}\) and transition probabilities \(\Gamma_{[a, b, k, l]}\). The stochastic_param roughly encodes the causal relationship between, on the one hand, \(U^{\text{group}}_{[g]}\), \(U^{\text{ind}}_{[a]}\), \(U^{\text{ind}}_{[b]}\), and \(U^{\text{dyad}}_{[a, b]}\), and, on the other hand, \(\theta_{[a, b, k]}\) and \(\Gamma_{[a, b, k, l]}\), shown at the top-left of our plated causal diagram (Figure 4).
We run scatter and density to display the distributions of all stochastic parameters and their pairwise association. (We do not display the resulting figure in the Quarto notebooks, but it can be found in the supplementary materials of the associated paper.)
Now that we have the parameters \(\theta\) and \(\Gamma\) for each dyad, we can specify how these parameters give rise to social interactions. The ctmc function produces the latent sequence of sojourns \(i\) for the dyads under study. The function takes as input: (i) the number of individuals per group; (ii) the amount of time the function should run for per group; (iii) params, generated by stochastic_param. As output, it produces a list, where the first level is group, and the second level is a data frame containing all of the the sojourns \(i\) for a given dyad. This function corresponds to Algorithm 1, in the supplemantary materials of the associated manuscript.
ctmc <-function(N_ind_vec =c(3, 3, 4),N_groups =length(N_ind_vec),N_dyad_vec = ((N_ind_vec * N_ind_vec) - N_ind_vec) /2, params,end_sim =c(1e3, 1e3, 1e3), ...) {# First level of df is group df <-list()for (group in1:N_groups) {# Second level of df is dyad df[[group]] <-list() N_dyad <- N_dyad_vec[group]# Next, we draw sojourns, each of which is made of two rows# one for its start, one for its endfor (dyad in1:N_dyad) {## 1. First sojourn (state 1) (df[[group]][[dyad]] <-tibble(id =c(1, 1), # observation idx =rexp(1, 1/ params[[group]]$theta[dyad, 1]) %>%rep(2), # holding timet =c(0, x[1]), # time stampss =c(1, 1), # statez =c("start", "end"), # modifierc =0# censoring level ))## 2. Run the CTMCrepeat { last_state <- df[[group]][[dyad]]$s %>%last()# Draw new state from a multinomial distribution new_state <-rmultinom(1, 1, params[[group]]$Gamma[, , dyad][last_state,]) %>%which.max()# Next sojourn df[[group]][[dyad]] <- df[[group]][[dyad]] %>%add_row(id = df[[group]][[dyad]]$id %>%last() +1,x =rexp(1, 1/ params[[group]]$theta[dyad, new_state]),t = df[[group]][[dyad]]$t %>%last(),s = new_state,z ="start",c =0 ) df[[group]][[dyad]] <- df[[group]][[dyad]] %>%add_row(id = df[[group]][[dyad]]$id %>%last(),x = df[[group]][[dyad]]$x %>%last(),t = (df[[group]][[dyad]]$t %>%last()) +last(df[[group]][[dyad]]$x),s = df[[group]][[dyad]]$s %>%last(),z ="end",c =0 )# break if we're after the maximal timeif (tail(df[[group]][[dyad]]$t, 1) >= end_sim[group]) {break } # end if } # end break } # end for dyad } # end for group df %>%return()} # end ctmc
For example, seq_states_ex[[1]][[2]], printed below, represents the sequence of sojourns that take place for the second dyad in the first group, up to a certain time limit \(t = 5000\) time units.
[[1]]
[[1]][[1]]
# A tibble: 10 × 6
id x t s z c
<dbl> <dbl> <dbl> <dbl> <chr> <dbl>
1 1 3041. 0 1 start 0
2 1 3041. 3041. 1 end 0
3 2 0.848 3041. 2 start 0
4 2 0.848 3041. 2 end 0
5 3 0.539 3041. 3 start 0
6 3 0.539 3042. 3 end 0
7 4 0.752 3042. 2 start 0
8 4 0.752 3043. 2 end 0
9 5 5063. 3043. 1 start 0
10 5 5063. 8105. 1 end 0
[[1]][[2]]
# A tibble: 18 × 6
id x t s z c
<dbl> <dbl> <dbl> <dbl> <chr> <dbl>
1 1 2490. 0 1 start 0
2 1 2490. 2490. 1 end 0
3 2 8.30 2490. 4 start 0
4 2 8.30 2498. 4 end 0
5 3 1961. 2498. 1 start 0
6 3 1961. 4459. 1 end 0
7 4 0.367 4459. 2 start 0
8 4 0.367 4460. 2 end 0
9 5 88.3 4460. 1 start 0
10 5 88.3 4548. 1 end 0
11 6 0.536 4548. 4 start 0
12 6 0.536 4549. 4 end 0
13 7 1.35 4549. 2 start 0
14 7 1.35 4550. 2 end 0
15 8 2.47 4550. 3 start 0
16 8 2.47 4552. 3 end 0
17 9 1261. 4552. 1 start 0
18 9 1261. 5813. 1 end 0
[[1]][[3]]
# A tibble: 16 × 6
id x t s z c
<dbl> <dbl> <dbl> <dbl> <chr> <dbl>
1 1 1156. 0 1 start 0
2 1 1156. 1156. 1 end 0
3 2 0.490 1156. 3 start 0
4 2 0.490 1157. 3 end 0
5 3 3484. 1157. 1 start 0
6 3 3484. 4641. 1 end 0
7 4 7.53 4641. 4 start 0
8 4 7.53 4648. 4 end 0
9 5 194. 4648. 1 start 0
10 5 194. 4842. 1 end 0
11 6 0.221 4842. 3 start 0
12 6 0.221 4842. 3 end 0
13 7 1.84 4842. 2 start 0
14 7 1.84 4844. 2 end 0
15 8 3282. 4844. 1 start 0
16 8 3282. 8126. 1 end 0
[[2]]
[[2]][[1]]
# A tibble: 4 × 6
id x t s z c
<dbl> <dbl> <dbl> <dbl> <chr> <dbl>
1 1 5369. 0 1 start 0
2 1 5369. 5369. 1 end 0
3 2 7.64 5369. 4 start 0
4 2 7.64 5377. 4 end 0
[[2]][[2]]
# A tibble: 6 × 6
id x t s z c
<dbl> <dbl> <dbl> <dbl> <chr> <dbl>
1 1 2013. 0 1 start 0
2 1 2013. 2013. 1 end 0
3 2 0.598 2013. 4 start 0
4 2 0.598 2013. 4 end 0
5 3 10201. 2013. 1 start 0
6 3 10201. 12214. 1 end 0
[[2]][[3]]
# A tibble: 6 × 6
id x t s z c
<dbl> <dbl> <dbl> <dbl> <chr> <dbl>
1 1 3786. 0 1 start 0
2 1 3786. 3786. 1 end 0
3 2 1.81 3786. 2 start 0
4 2 1.81 3788. 2 end 0
5 3 9421. 3788. 1 start 0
6 3 9421. 13209. 1 end 0
[[2]][[4]]
# A tibble: 32 × 6
id x t s z c
<dbl> <dbl> <dbl> <dbl> <chr> <dbl>
1 1 340. 0 1 start 0
2 1 340. 340. 1 end 0
3 2 5.15 340. 4 start 0
4 2 5.15 345. 4 end 0
5 3 715. 345. 1 start 0
6 3 715. 1060. 1 end 0
7 4 8.93 1060. 3 start 0
8 4 8.93 1069. 3 end 0
9 5 7.08 1069. 4 start 0
10 5 7.08 1076. 4 end 0
# ℹ 22 more rows
[[2]][[5]]
# A tibble: 4 × 6
id x t s z c
<dbl> <dbl> <dbl> <dbl> <chr> <dbl>
1 1 9557. 0 1 start 0
2 1 9557. 9557. 1 end 0
3 2 4.97 9557. 3 start 0
4 2 4.97 9562. 3 end 0
[[2]][[6]]
# A tibble: 4 × 6
id x t s z c
<dbl> <dbl> <dbl> <dbl> <chr> <dbl>
1 1 7298. 0 1 start 0
2 1 7298. 7298. 1 end 0
3 2 5.42 7298. 4 start 0
4 2 5.42 7303. 4 end 0
1.4d_ff
Now that we can generate the true dyadic behavioural sojourns with ctmc, we can introduce measurement. d_ff does this by creating a set of dyadic focal follow, i.e. it cuts a 40 minutes windows from the latent sequence of states generated with ctmc for each dyad under consideration. It takes as inputs: (i) The output from ctmc: a list of data frames, where each data frame encodes the true sequences of sojourns \(i\). (ii) Two time stamps: the starting time and the end time for a 40 min focal follow. As outputs, the function returns a list of data frames, where each data frame contains the sojourns \(j\) that are basically the observations \(i\) from input (i) in between the time stamps given by input (ii). d_ff corresponds to Algorithm 2, in the supplemantary materials of the associated manuscript.
Note that in the next two functions, the censoring level c is related to the variable \(C\) described in the paper, but somewhat deviates from it.
d_ff <-function(df, start_focal, end_focal){for (dyad in1:length(df)){# Censor the start of the dyadic focal follow # i.e. we cut anything before start_focal df[[dyad]] <- df[[dyad]] %>%add_row(t = start_focal, c =1) %>%# c = 1 means we don't know the star timeadd_row(t = end_focal, c =2) %>%#c = 2 means we don't know the end timearrange(t) %>%mutate(id =case_when( c ==1~lag(id), c ==2~lead(id),TRUE~ id ),s =case_when( c ==1~lag(s), c ==2~lead(s),TRUE~ s ),z =case_when( c ==1~"start", c ==2~"end",TRUE~ z ),x =case_when( c ==1~lead(t) - start_focal, c ==2~ end_focal -lag(t),TRUE~ x ) ) %>%group_by(id) %>%mutate(x =case_when(sum(c) ==0~ x, sum(c) %in%c(1, 2) ~ x[2],sum(c) ==3~ end_focal - start_focal ) ) %>%ungroup() %>%filter(t >= start_focal & t <= end_focal) } # end for ireturn(df)} # end d_ff
For example:
# Before cutting seq_states_ex[[2]][5:6]
[[1]]
# A tibble: 4 × 6
id x t s z c
<dbl> <dbl> <dbl> <dbl> <chr> <dbl>
1 1 9557. 0 1 start 0
2 1 9557. 9557. 1 end 0
3 2 4.97 9557. 3 start 0
4 2 4.97 9562. 3 end 0
[[2]]
# A tibble: 4 × 6
id x t s z c
<dbl> <dbl> <dbl> <dbl> <chr> <dbl>
1 1 7298. 0 1 start 0
2 1 7298. 7298. 1 end 0
3 2 5.42 7298. 4 start 0
4 2 5.42 7303. 4 end 0
# After cutting seq_states_ex[[2]][5:6] %>%d_ff(start_focal =820, end_focal =860)
[[1]]
# A tibble: 2 × 6
id x t s z c
<dbl> <dbl> <dbl> <dbl> <chr> <dbl>
1 1 40 820 1 start 1
2 1 40 860 1 end 2
[[2]]
# A tibble: 2 × 6
id x t s z c
<dbl> <dbl> <dbl> <dbl> <chr> <dbl>
1 1 40 820 1 start 1
2 1 40 860 1 end 2
1.5sim
The sim function takes, as input, a set of stochastic parameters generated with stochastic_param, as well as the number of focal follows, burnin period, and focal duration. It produces, as output, behavioural observations that can be fitted to a Stan model.
sim <-function(N_ind_vec,N_group =length(N_ind_vec),N_dyad_vec = ((N_ind_vec * N_ind_vec) - N_ind_vec) /2, parameters,n_focals_vec = N_ind_vec *10,burnin =1e3,duration_focal =40,end_sim_vec = burnin + n_focals_vec * duration_focal) {## 1. Define a number of objects dyads <-list() # list of dyads and individuals per group sim_data <-list() # simulated observations (sojourns in states per dyad) d <-list() # summary of sim_data## 2. We run ctmc for each group state_seq <-ctmc(N_ind_vec = N_ind_vec,params = parameters,end_sim = end_sim_vec)for (group in1:N_group) {## 3. Define individuals & dyads# We give individuals unique name across groups# e.g. 1:10 (group 1), 11:20 (group 2), etc. bl_ind <-c(0, head(N_ind_vec, -1))[1:group] %>%sum() bl_dyad <-c(0, head(N_dyad_vec, -1))[1:group] %>%sum() N_ind <- N_ind_vec[group] dyads[[group]] <-tibble(grp = group,ind_a = bl_ind +t(combn(N_ind, 2))[, 1],ind_b = bl_ind +t(combn(N_ind, 2))[, 2],dyad_id = bl_dyad +c(1:N_dyad_vec[group]) )## 4. Define sequence of focal follows for each group# 4.1 Sequence of individuals ids to focal by group# We make sure to focal each individual at least once focal_ind <- bl_ind +c(1:N_ind)# For the remaining focal-follows, we draw individuals at random focal_ind[(N_ind +1):n_focals_vec[group]] <-sample(x = bl_ind +1:N_ind,size = n_focals_vec[group] - N_ind,replace =TRUE)# 4.2 We mark the start time of the focal-follows time_stamps <-seq(from = burnin, to = end_sim_vec[[group]], by = duration_focal) dyads_to_sample <-list() # dyads to sample for each focal follow sim_data[[group]] <-list() # observed sequence of states d[[group]] <-list()# for each individual focal followfor (focal in1:n_focals_vec[group]) {# sim_data is indexed by (i) group; (ii) ind. focal follow; (iii) dyad. sim_data[[group]][[focal]] <-list() d[[group]][[focal]] <-list()## 5. We conduct the focal follows per se# Which *dyads* are we practically following for each ind. focal-follow dyads_to_sample[[focal]] <- dyads[[group]] %>%filter(ind_a == focal_ind[focal] | ind_b == focal_ind[focal]) %>%pull(dyad_id)# We select the latent state sequence of interest sim_data[[group]][[focal]] <- state_seq[[group]][dyads_to_sample[[focal]] - bl_dyad]# We cut the windows corresponding to focal-follows sim_data[[group]][[focal]] <-d_ff(sim_data[[group]][[focal]], time_stamps[focal], time_stamps[focal] + duration_focal)# for each dyadfor (dyad in1:(N_ind -1)) {## add name of focal individual and focal dyad for each sim_data[[group]][[focal]][[dyad]] <- sim_data[[group]][[focal]][[dyad]] %>%mutate(focal_id = focal_ind[focal],focal_dyad = dyads_to_sample[[focal]][dyad])## 6. Summarise each sojourn by a few variables d[[group]][[focal]][[dyad]] <- sim_data[[group]][[focal]][[dyad]] %>%group_by(id) %>%# for each state sojourn idsummarize(grp = group, # grp idx =first(x), # holding times =first(s), # statec =ifelse(sum(c) <2, 0, 1),# censoring = {0, 1}dyad =first(focal_dyad), # id of dyadfocal_id =first(focal_id), # id of focal individualprot_id = focal # id of the protocol ) %>%select(-id) } # end dyad } # end for focal d[[group]] <- d[[group]] %>%bind_rows() } # end for group## 7. We collapse the sojourn data d <- d %>%bind_rows() %>%group_by(grp, focal_id, dyad, prot_id) %>%# add a unique identifier per dyadic focal followmutate(prot_id =cur_group_id()) %>%ungroup() %>%arrange(prot_id)## 8. Transitions obs_tr <-tibble(grp =head(d$grp, -1),dyad_id =head(d$dyad, -1),prot_id_from =head(d$prot_id, -1),prot_id_to =tail(d$prot_id, -1),s_from =head(d$s, -1),s_to =tail(d$s, -1) ) %>%filter(prot_id_from == prot_id_to) %>%select(-c(prot_id_from, prot_id_to))## 9. Exposure and "events" in each state exp_events <- d %>%mutate(observed =ifelse(c ==0, 1, 0)) %>%group_by(dyad, s) %>%# Amount of time spent in each statesummarise(exposure =sum(x),# how many times have we observed the end?events =sum(observed)) %>%right_join(expand.grid(dyad =1:sum(N_dyad_vec),s =1:4)) %>%arrange(dyad, s) %>%mutate_at(c("exposure", "events"), ~replace_na(., 0))# Stan-friendly format list <-list()# General indices list$N_dyad <-sum(N_dyad_vec) list$N_ind <-sum(N_ind_vec) list$N_group <- N_group list$N_ind_per_group <- N_ind_vec list$N_dyad_per_group <- N_dyad_vec # Individuals and dyads (dataset "A") list$A_grp <-bind_rows(dyads)$grp list$A_ind_a <-bind_rows(dyads)$ind_a list$A_ind_b <-bind_rows(dyads)$ind_b list$A_dyad <-bind_rows(dyads)$dyad_id# Observations j: sequence of states (dataset "B") list$J <-bind_rows(d) %>%nrow() list$B_grp <- d$grp list$B_x <- d$x list$B_s <- d$s list$B_c <- d$c list$B_dyad <- d$dyad# Observed transitions (dataset "C") list$N_trans <-nrow(obs_tr) list$C_grp <- obs_tr$grp list$C_dyad <- obs_tr$dyad_id list$C_s_from <- obs_tr$s_from list$C_s_to <- obs_tr$s_to# Observed exposure time and count of events in each state (dataset "D") list$N_row_D <-nrow(exp_events) list$D_dyad <- exp_events$dyad list$D_s <- exp_events$s list$D_exposure <- exp_events$exposure list$D_events <- exp_events$events list %>%return() } # end sim
1.6 Run sim
We run sim for four groups of ten individuals, conducting \(1500\) focal follows for each group.
param_1 <-list() # stochastic parametersd_1 <-list() # simulated dataset.seed(247)for (i in1:10) { param_1[[i]] <-stochastic_param(N_ind_vec =c(10, 10, 10, 10)) d_1[[i]] <-sim(N_ind_vec =c(10, 10, 10, 10),n_focals_vec =c(1500, 1500, 1500, 1500),parameters = param_1[[i]])}# We save the data param_1 %>%saveRDS("./sim_data/01/param_1.rds") d_1 %>%saveRDS("./sim_data/01/d_1.rds")
2 Observed distributions
Import the generate data:
# Import data d_1 <-readRDS("./sim_data/01/d_1.rds")
We plot the number of observations \(j\) per dyad \((a, b)\) and per state \(k\).
For instance, consider the panel on the left-hand side (S = 1). We observe that a high number of dyads have about \(\simeq 300\) observations \(j\) in state \(k = 1\). In contrast, very few dyads have more than \(\simeq 340\) samples in state \(k = 1\).
Now, let us plot the number of observed transitions per transition types \(k \to l\) and per dyad \((a, b)\):
# Stan model (m_1 <-cmdstan_model("./stan_models/core_model.stan"))
data{
// Indices
int<lower = 1> N_dyad;
int<lower = 1> N_ind;
int<lower = 1> N_group;
// Dataset "A": all dyads, and corresponding grp and ind
array[N_dyad] int A_dyad; // Dyad id {1, ..., N_dyad}
array[N_dyad] int A_grp; // Group id {1, ..., N_group}
array[N_dyad] int A_ind_a; // ind a per dyad number
array[N_dyad] int A_ind_b; // ind b per dyad number
// Dataset "B": observed states
int<lower = 1> J; // Number of observed states
array[J] int<lower = 1, upper = N_group> B_grp; // Group
vector<lower = 0>[J] B_x; // Holding time
array[J] int<lower = 1, upper = 4> B_s; // State {1, 2, 3, 4}
array[J] int<lower = 0, upper = 1> B_c; // censoring {0, 1}
array[J] int<lower = 1, upper = N_dyad> B_dyad; // dyad {1, ..., N_dyad}
// Dataset "C": observed transitions
int<lower = 1> N_trans; // Number of observed transitions
array[N_trans] int<lower = 1, upper = N_group> C_grp; // Group
array[N_trans] int<lower = 1, upper = N_dyad> C_dyad; // dyad {1, ..., N_dyad}
array[N_trans] int<lower = 1, upper = 4> C_s_from; // state j
array[N_trans] int<lower = 1, upper = 4> C_s_to; // state j + 1
}
transformed data {
// Count of transitions. One matrix per dyad
array[N_dyad, 4, 4] int<lower = 0> transition_counts;
// Initialise the transition matrices with zeros
for (dyad in 1:N_dyad) {
for (k in 1:4) {
for (l in 1:4) {
transition_counts[dyad, k, l] = 0;
}
}
}
// Fill the transition matrices with transitions from the input data
for (tr in 1:N_trans) {
transition_counts[C_dyad[tr], C_s_from[tr], C_s_to[tr]] += 1;
}
}
parameters {
// Intercepts per group
vector[N_group] alpha_1;
vector[N_group] alpha_2;
vector[N_group] alpha_gr;
vector[N_group] alpha_1_gr;
vector[N_group] alpha_2_gr;
vector[N_group] alpha_gr_1;
vector[N_group] alpha_gr_2;
// Individual varying effects
matrix[12, N_ind] z_phi; // z-score of phi (wide format: 12 R, N_ind C)
vector<lower = 0>[12] sigma_phi; // SD of phi
cholesky_factor_corr[12] L_phi; // Cholesky factor of corr. matrix phi
// Dyadic varying effects
matrix[12, N_dyad] z_tau; // z-score of tau (wide format: 12 R, N_dyad C)
vector<lower = 0>[7] sigma_tau; // SD of tau
cholesky_factor_corr[12] L_tau; // Cholesky factor of corr. matrix for tau
}
transformed parameters{
// 1. individual varying effects
matrix[N_ind, 12] phi_raw; // vertical: N_ind R, 12 C)
phi_raw = (diag_pre_multiply(sigma_phi, L_phi) * z_phi)';
vector[N_ind] phi_1 = phi_raw[, 1];
vector[N_ind] phi_2 = phi_raw[, 2];
vector[N_ind] phi_give_gr = phi_raw[, 3];
vector[N_ind] phi_rec_gr = phi_raw[, 4];
vector[N_ind] phi_1_give_gr = phi_raw[, 5];
vector[N_ind] phi_1_rec_gr = phi_raw[, 6];
vector[N_ind] phi_2_give_gr = phi_raw[, 7];
vector[N_ind] phi_2_rec_gr = phi_raw[, 8];
vector[N_ind] phi_give_gr_1 = phi_raw[, 9];
vector[N_ind] phi_give_gr_2 = phi_raw[, 10];
vector[N_ind] phi_rec_gr_1 = phi_raw[, 11];
vector[N_ind] phi_rec_gr_2 = phi_raw[, 12];
// 2. dyad varying effects
matrix[N_dyad, 12] tau_raw; // vertical: N_ind R, 12 C
tau_raw = (diag_pre_multiply([sigma_tau[1], sigma_tau[2],
sigma_tau[3], sigma_tau[3],
sigma_tau[4], sigma_tau[4],
sigma_tau[5], sigma_tau[5],
sigma_tau[6], sigma_tau[7],
sigma_tau[6], sigma_tau[7]], L_tau) * z_tau)';
vector[N_dyad] tau_1 = tau_raw[, 1];
vector[N_dyad] tau_2 = tau_raw[, 2];
vector[N_dyad] tau_gr_ab = tau_raw[, 3];
vector[N_dyad] tau_gr_ba = tau_raw[, 4];
vector[N_dyad] tau_1_gr_ab = tau_raw[, 5];
vector[N_dyad] tau_1_gr_ba = tau_raw[, 6];
vector[N_dyad] tau_2_gr_ab = tau_raw[, 7];
vector[N_dyad] tau_2_gr_ba = tau_raw[, 8];
vector[N_dyad] tau_gr_1_ab = tau_raw[, 9];
vector[N_dyad] tau_gr_2_ab = tau_raw[, 10];
vector[N_dyad] tau_gr_1_ba = tau_raw[, 11];
vector[N_dyad] tau_gr_2_ba = tau_raw[, 12];
// We declare a bunch of varibales that we define further below
array[N_dyad, 4] real theta;
array[N_dyad, 4, 4] real Psi;
array[N_dyad, 4, 4] real Gamma;
// 3. theta
for (dyad in 1:N_dyad) {
theta[dyad, 1] = exp(alpha_1[A_grp[dyad]] +
phi_1[A_ind_a[dyad]] +
phi_1[A_ind_b[dyad]] +
tau_1[A_dyad[dyad]]);
theta[dyad, 2] = exp(alpha_2[A_grp[dyad]] +
phi_2[A_ind_a[dyad]] +
phi_2[A_ind_b[dyad]] +
tau_2[A_dyad[dyad]]);
theta[dyad, 3] = exp(alpha_gr[A_grp[dyad]] +
phi_give_gr[A_ind_a[dyad]] +
phi_rec_gr[A_ind_b[dyad]] +
tau_gr_ab[A_dyad[dyad]]);
theta[dyad, 4] = exp(alpha_gr[A_grp[dyad]] +
phi_give_gr[A_ind_b[dyad]] +
phi_rec_gr[A_ind_a[dyad]] +
tau_gr_ba[A_dyad[dyad]]);
// 4. Psi
Psi[dyad, 1, 2] = 0.0;
Psi[dyad, 1, 3] = alpha_1_gr[A_grp[dyad]] + phi_1_give_gr[A_ind_a[dyad]] +
phi_1_rec_gr[A_ind_b[dyad]] + tau_1_gr_ab[A_dyad[dyad]];
Psi[dyad, 1, 4] = alpha_1_gr[A_grp[dyad]] + phi_1_give_gr[A_ind_b[dyad]] +
phi_1_rec_gr[A_ind_a[dyad]] + tau_1_gr_ba[A_dyad[dyad]];
Psi[dyad, 2, 1] = 0.0;
Psi[dyad, 2, 3] = alpha_2_gr[A_grp[dyad]] + phi_2_give_gr[A_ind_a[dyad]] +
phi_2_rec_gr[A_ind_b[dyad]] + tau_2_gr_ab[A_dyad[dyad]];
Psi[dyad, 2, 4] = alpha_2_gr[A_grp[dyad]] + phi_2_give_gr[A_ind_b[dyad]] +
phi_2_rec_gr[A_ind_a[dyad]] + tau_2_gr_ba[A_dyad[dyad]];
Psi[dyad, 3, 1] = alpha_gr_1[A_grp[dyad]] + phi_give_gr_1[A_ind_a[dyad]] +
phi_rec_gr_1[A_ind_b[dyad]] + tau_gr_1_ab[A_dyad[dyad]];
Psi[dyad, 3, 2] = alpha_gr_2[A_grp[dyad]] + phi_give_gr_2[A_ind_a[dyad]] +
phi_rec_gr_2[A_ind_b[dyad]] + tau_gr_2_ab[A_dyad[dyad]];
Psi[dyad, 3, 4] = 0.0;
Psi[dyad, 4, 1] = alpha_gr_1[A_grp[dyad]] + phi_give_gr_1[A_ind_b[dyad]] +
phi_rec_gr_1[A_ind_a[dyad]] + tau_gr_1_ba[A_dyad[dyad]];
Psi[dyad, 4, 2] = alpha_gr_2[A_grp[dyad]] + phi_give_gr_2[A_ind_b[dyad]] +
phi_rec_gr_2[A_ind_a[dyad]] + tau_gr_2_ba[A_dyad[dyad]];
Psi[dyad, 4, 3] = 0.0;
// 5. Gamma
for (k in 1:4) {
real row_sum = 0.0;
// Compute the denominator (sum of exponents in row k, excluding diagonal)
for (m in 1:4) {
if (m != k) {
row_sum += exp(Psi[dyad, k, m]);
} // end if
} // end for m
// Now compute Gamma for each element in row k
for (l in 1:4) {
if (k == l) {
Gamma[dyad, k, l] = 0.0; // Diagonal elements are 0
} else {
Gamma[dyad, k, l] = exp(Psi[dyad, k, l]) / row_sum;
} // end if
} // end for l
} // end for k
} // end for dyad
}
model{
// 1. Likelihood for holding times
for (j in 1:J){
// not censored
if (B_c[j] == 0){
target += exponential_lpdf(B_x[j] | 1 / theta[B_dyad[j], B_s[j]]);
} // end non-censored obs
// censored
if (B_c[j] == 1){
target += exponential_lccdf(B_x[j] | 1 / theta[B_dyad[j], B_s[j]]);
} // end censored obs
} // end for j
// 2. Likelihood for observed transitions
for (tr in 1:N_trans) {
target += categorical_lpmf(C_s_to[tr] |
to_vector(Gamma[C_dyad[tr], C_s_from[tr]]));
} // end for tr
// 3. Hyper-priors
target += normal_lpdf(alpha_1 | 8, 2);
target += normal_lpdf(alpha_2 | 1.5, 1);
target += normal_lpdf(alpha_gr | 1.5, 1);
target += normal_lpdf(alpha_1_gr | 0, 1);
target += normal_lpdf(alpha_2_gr | 0, 1);
target += normal_lpdf(alpha_gr_1 | 0, 1);
target += normal_lpdf(alpha_gr_2 | 0, 1);
target += normal_lpdf(to_vector(z_phi) | 0, 1);
target += normal_lpdf(to_vector(z_tau) | 0, 1);
target += exponential_lpdf(sigma_phi | 1);
target += exponential_lpdf(sigma_tau | 1);
target += lkj_corr_cholesky_lpdf(L_phi | 4);
target += lkj_corr_cholesky_lpdf(L_tau | 4);
} // end model block
generated quantities{
// Corr. matrix from cholesky factors
// Varying individual effects
matrix[12, 12] c_ind;
c_ind = multiply_lower_tri_self_transpose(L_phi);
// Varying dyadic effects
matrix[12, 12] c_dyad;
c_dyad = multiply_lower_tri_self_transpose(L_tau);
// Avg. holding time per state (across dyads)
vector[3] avg_theta;
avg_theta[1] = mean(theta[, 1]);
avg_theta[2] = mean(theta[, 2]);
avg_theta[3] = mean(append_row(
to_vector(theta[, 3]), to_vector(theta[, 4])));
// Avg. transition probabilities (across dyads)
matrix[3, 3] avg_Gamma;
avg_Gamma[1, 2] = mean(Gamma[, 1, 2]);
avg_Gamma[1, 3] = mean(Gamma[, 1, 3]) + mean(Gamma[, 1, 4]);
avg_Gamma[2, 1] = mean(Gamma[, 2, 1]);
avg_Gamma[2, 3] = mean(Gamma[, 2, 3]) + mean(Gamma[, 2, 4]);
avg_Gamma[3, 1] = mean(append_row(
to_vector(Gamma[, 3, 1]), to_vector(Gamma[, 4, 1])));
avg_Gamma[3, 2] = mean(append_row(
to_vector(Gamma[, 3, 2]), to_vector(Gamma[, 4, 2])));
avg_Gamma[3, 3] = mean(append_row(
to_vector(Gamma[, 3, 4]), to_vector(Gamma[, 4, 3])));
}
3.1 MCMC
We run the Stan models, obtaining MCMC samples of the posterior distribution.
# We define lists post_1 <-list()# Import data d_1 <-readRDS("./sim_data/01/d_1.rds")# Model without transformed parametersfor(i in1:10){ post_1[[i]] <- m_1$sample(data = d_1[[i]],iter_warmup =1e3,iter_sampling =1e3,chains =4,parallel_chains =4 )if(i ==1){diag_1 <- post_1[[i]]$summary()}# Extract posterior samples post_1[[i]] <- post_1[[i]] %>%tidy_draws() %>%mutate(iter = i) }# We reformat the posterior and save it post_1 %>%bind_rows() %>%saveRDS("./fitted_models/01.1/post_1.rds") diag_1 %>%saveRDS("./fitted_models/01.1/diag_1.rds")
We import the posterior samples (post_1), posterior diagnostics (diag_1), and stochastic parameters (param_1).
# Import posterior draws and parameter values post_1 <-readRDS("./fitted_models/01.1/post_1.rds") param_1 <-readRDS("./sim_data/01/param_1.rds") diag_1 <-readRDS("./fitted_models/01.1/diag_1.rds")
3.2 Theta and gamma
In this section, we plot the estimated varying effects \(\theta_{[a, b, k]}\) and \(\Gamma_{[a, b, k, l]}\) and compare to the target values from the generative model.
We define the dyads under scrutiny:
# Iteration of the generative modeliteration <-1# Sample 10 dyadsset.seed(3666)dyads <-sample(1:180, 10, replace =FALSE)# These dyads are sampled across all four groupsgrp_nbs <-1:4# Rename and import relevant objectspost <- post_1data <-readRDS("./sim_data/01/d_1.rds")parameters <- param_1
Below, we average \(\theta\) and \(\Gamma\) over the three groups for a number of iterations, then average over all the previous iterations. By doing so, we compute the population averages for holding times and transition probabilities.
For each posterior, and for each observed sojourn in a given state, we draw a posterior prediction. We then only keep sojourns smaller than 40 minutes.
# We bin the observed holding times for plottingpp_x_tbl <- post_pred_x %>%bind_rows() %>%# Binning sequencemutate(bin =sapply(x, function(val)seq(2.5, 37.5, by =5)[which.min(abs(seq(2.5, 37.5, by =5) - val))])) %>%group_by(s, sample, bin) %>%summarise(count =n()) %>%group_by(s, sample) %>%# Fill in the zeroscomplete(bin =seq(2.5, 37.5, by =5), fill =list(count =0)) %>%ungroup()p1 <- obs_ht %>%filter(s ==1& c ==0) %>%hist_plot(colbin ="#c9c7bd",posterior =1,post_data = pp_x_tbl %>%filter(s ==1),alpha_point =0.6,alpha_line =0.4,alpha_hist =0.4,col_post ="#c9c7bd" )p2 <- obs_ht %>%filter(s ==2& c %in%c(0, 1)) %>%hist_plot(colbin ="#E9C4C1",posterior =1,post_data = pp_x_tbl %>%filter(s ==2),alpha_point =0.6,alpha_line =0.4,alpha_hist =0.4,col_post ="#d1a7a3" )p3 <- obs_ht %>%filter(s ==3& c %in%c(0, 1)) %>%hist_plot(colbin ="#ce8da6",posterior =1,post_data = pp_x_tbl %>%filter(s ==3),alpha_point =0.6,alpha_line =0.4,alpha_hist =0.4,col_post ="#a37184" )p4 <- obs_ht %>%filter(s ==4& c %in%c(0, 1)) %>%hist_plot(colbin ="#996282",posterior =1,post_data = pp_x_tbl %>%filter(s ==4),alpha_point =0.6,alpha_line =0.4,alpha_hist =0.4,col_post ="#7a4e68" )
In this section, we show the prior distribution of parameters \(\theta\) and \(\Gamma\), and of the observables \(X\).
m_prior differs from m_1 in two ways. First, \(\theta\) and \(\Gamma\) are defined as transformed parameters, so we get posterior samples for them. Second, there is no likelihood function, so the posterior distribution is nothing more than the prior distribution.
# Stan model (m_prior <-cmdstan_model("./stan_models/prior_core_model.stan"))
data{
// Indices
int<lower = 1> N_dyad;
int<lower = 1> N_ind;
int<lower = 1> N_group;
// Dataset "A": all dyads, and corresponding grp and ind
array[N_dyad] int A_dyad; // Dyad id {1, ..., N_dyad}
array[N_dyad] int A_grp; // Group id {1, ..., N_group}
array[N_dyad] int A_ind_a; // ind a per dyad number
array[N_dyad] int A_ind_b; // ind b per dyad number
// Dataset "B": observed states
int<lower = 1> J; // Number of observed states
array[J] int<lower = 1, upper = N_group> B_grp; // Group
vector<lower = 0>[J] B_x; // Holding time
array[J] int<lower = 1, upper = 4> B_s; // State {1, 2, 3, 4}
array[J] int<lower = 0, upper = 1> B_c; // censoring {0, 1}
array[J] int<lower = 1, upper = N_dyad> B_dyad; // dyad {1, ..., N_dyad}
// Dataset "C": observed transitions
int<lower = 1> N_trans; // Number of observed transitions
array[N_trans] int<lower = 1, upper = N_group> C_grp; // Group
array[N_trans] int<lower = 1, upper = N_dyad> C_dyad; // dyad {1, ..., N_dyad}
array[N_trans] int<lower = 1, upper = 4> C_s_from; // state j
array[N_trans] int<lower = 1, upper = 4> C_s_to; // state j + 1
}
transformed data {
// Count of transitions. One matrix per dyad
array[N_dyad, 4, 4] int<lower = 0> transition_counts;
// Initialise the transition matrices with zeros
for (dyad in 1:N_dyad) {
for (k in 1:4) {
for (l in 1:4) {
transition_counts[dyad, k, l] = 0;
}
}
}
// Fill the transition matrices with transitions from the input data
for (tr in 1:N_trans) {
transition_counts[C_dyad[tr], C_s_from[tr], C_s_to[tr]] += 1;
}
}
parameters {
// Intercepts per group
vector[N_group] alpha_1;
vector[N_group] alpha_2;
vector[N_group] alpha_gr;
vector[N_group] alpha_1_gr;
vector[N_group] alpha_2_gr;
vector[N_group] alpha_gr_1;
vector[N_group] alpha_gr_2;
// Individual varying effects
matrix[12, N_ind] z_phi; // z-score of phi (wide format: 12 R, N_ind C)
vector<lower = 0>[12] sigma_phi; // SD of phi
cholesky_factor_corr[12] L_phi; // Cholesky factor of corr. matrix phi
// Dyadic varying effects
matrix[12, N_dyad] z_tau; // z-score of tau (wide format: 12 R, N_dyad C)
vector<lower = 0>[7] sigma_tau; // SD of tau
cholesky_factor_corr[12] L_tau; // Cholesky factor of corr. matrix for tau
}
transformed parameters{
// 1. individual varying effects
matrix[N_ind, 12] phi_raw; // vertical: N_ind R, 12 C)
phi_raw = (diag_pre_multiply(sigma_phi, L_phi) * z_phi)';
vector[N_ind] phi_1 = phi_raw[, 1];
vector[N_ind] phi_2 = phi_raw[, 2];
vector[N_ind] phi_give_gr = phi_raw[, 3];
vector[N_ind] phi_rec_gr = phi_raw[, 4];
vector[N_ind] phi_1_give_gr = phi_raw[, 5];
vector[N_ind] phi_1_rec_gr = phi_raw[, 6];
vector[N_ind] phi_2_give_gr = phi_raw[, 7];
vector[N_ind] phi_2_rec_gr = phi_raw[, 8];
vector[N_ind] phi_give_gr_1 = phi_raw[, 9];
vector[N_ind] phi_give_gr_2 = phi_raw[, 10];
vector[N_ind] phi_rec_gr_1 = phi_raw[, 11];
vector[N_ind] phi_rec_gr_2 = phi_raw[, 12];
// 2. dyad varying effects
matrix[N_dyad, 12] tau_raw; // vertical: N_ind R, 12 C
tau_raw = (diag_pre_multiply([sigma_tau[1], sigma_tau[2],
sigma_tau[3], sigma_tau[3],
sigma_tau[4], sigma_tau[4],
sigma_tau[5], sigma_tau[5],
sigma_tau[6], sigma_tau[7],
sigma_tau[6], sigma_tau[7]], L_tau) * z_tau)';
vector[N_dyad] tau_1 = tau_raw[, 1];
vector[N_dyad] tau_2 = tau_raw[, 2];
vector[N_dyad] tau_gr_ab = tau_raw[, 3];
vector[N_dyad] tau_gr_ba = tau_raw[, 4];
vector[N_dyad] tau_1_gr_ab = tau_raw[, 5];
vector[N_dyad] tau_1_gr_ba = tau_raw[, 6];
vector[N_dyad] tau_2_gr_ab = tau_raw[, 7];
vector[N_dyad] tau_2_gr_ba = tau_raw[, 8];
vector[N_dyad] tau_gr_1_ab = tau_raw[, 9];
vector[N_dyad] tau_gr_2_ab = tau_raw[, 10];
vector[N_dyad] tau_gr_1_ba = tau_raw[, 11];
vector[N_dyad] tau_gr_2_ba = tau_raw[, 12];
// We declare a bunch of varibales that we define further below
array[N_dyad, 4] real theta;
array[N_dyad, 4, 4] real Psi;
array[N_dyad, 4, 4] real Gamma;
// 3. theta
for (dyad in 1:N_dyad) {
theta[dyad, 1] = exp(alpha_1[A_grp[dyad]] +
phi_1[A_ind_a[dyad]] +
phi_1[A_ind_b[dyad]] +
tau_1[A_dyad[dyad]]);
theta[dyad, 2] = exp(alpha_2[A_grp[dyad]] +
phi_2[A_ind_a[dyad]] +
phi_2[A_ind_b[dyad]] +
tau_2[A_dyad[dyad]]);
theta[dyad, 3] = exp(alpha_gr[A_grp[dyad]] +
phi_give_gr[A_ind_a[dyad]] +
phi_rec_gr[A_ind_b[dyad]] +
tau_gr_ab[A_dyad[dyad]]);
theta[dyad, 4] = exp(alpha_gr[A_grp[dyad]] +
phi_give_gr[A_ind_b[dyad]] +
phi_rec_gr[A_ind_a[dyad]] +
tau_gr_ba[A_dyad[dyad]]);
// 4. Psi
Psi[dyad, 1, 2] = 0.0;
Psi[dyad, 1, 3] = alpha_1_gr[A_grp[dyad]] + phi_1_give_gr[A_ind_a[dyad]] +
phi_1_rec_gr[A_ind_b[dyad]] + tau_1_gr_ab[A_dyad[dyad]];
Psi[dyad, 1, 4] = alpha_1_gr[A_grp[dyad]] + phi_1_give_gr[A_ind_b[dyad]] +
phi_1_rec_gr[A_ind_a[dyad]] + tau_1_gr_ba[A_dyad[dyad]];
Psi[dyad, 2, 1] = 0.0;
Psi[dyad, 2, 3] = alpha_2_gr[A_grp[dyad]] + phi_2_give_gr[A_ind_a[dyad]] +
phi_2_rec_gr[A_ind_b[dyad]] + tau_2_gr_ab[A_dyad[dyad]];
Psi[dyad, 2, 4] = alpha_2_gr[A_grp[dyad]] + phi_2_give_gr[A_ind_b[dyad]] +
phi_2_rec_gr[A_ind_a[dyad]] + tau_2_gr_ba[A_dyad[dyad]];
Psi[dyad, 3, 1] = alpha_gr_1[A_grp[dyad]] + phi_give_gr_1[A_ind_a[dyad]] +
phi_rec_gr_1[A_ind_b[dyad]] + tau_gr_1_ab[A_dyad[dyad]];
Psi[dyad, 3, 2] = alpha_gr_2[A_grp[dyad]] + phi_give_gr_2[A_ind_a[dyad]] +
phi_rec_gr_2[A_ind_b[dyad]] + tau_gr_2_ab[A_dyad[dyad]];
Psi[dyad, 3, 4] = 0.0;
Psi[dyad, 4, 1] = alpha_gr_1[A_grp[dyad]] + phi_give_gr_1[A_ind_b[dyad]] +
phi_rec_gr_1[A_ind_a[dyad]] + tau_gr_1_ba[A_dyad[dyad]];
Psi[dyad, 4, 2] = alpha_gr_2[A_grp[dyad]] + phi_give_gr_2[A_ind_b[dyad]] +
phi_rec_gr_2[A_ind_a[dyad]] + tau_gr_2_ba[A_dyad[dyad]];
Psi[dyad, 4, 3] = 0.0;
// 5. Gamma
for (k in 1:4) {
real row_sum = 0.0;
// Compute the denominator (sum of exponents in row k, excluding diagonal)
for (m in 1:4) {
if (m != k) {
row_sum += exp(Psi[dyad, k, m]);
} // end if
} // end for m
// Now compute Gamma for each element in row k
for (l in 1:4) {
if (k == l) {
Gamma[dyad, k, l] = 0.0; // Diagonal elements are 0
} else {
Gamma[dyad, k, l] = exp(Psi[dyad, k, l]) / row_sum;
} // end if
} // end for l
} // end for k
} // end for dyad
}
model{
// Hyper-priors without likelihood
target += normal_lpdf(alpha_1 | 8, 2);
target += normal_lpdf(alpha_2 | 1.5, 1);
target += normal_lpdf(alpha_gr | 1.5, 1);
target += normal_lpdf(alpha_1_gr | 0, 1);
target += normal_lpdf(alpha_2_gr | 0, 1);
target += normal_lpdf(alpha_gr_1 | 0, 1);
target += normal_lpdf(alpha_gr_2 | 0, 1);
target += normal_lpdf(to_vector(z_phi) | 0, 1);
target += normal_lpdf(to_vector(z_tau) | 0, 1);
target += exponential_lpdf(sigma_phi | 1);
target += exponential_lpdf(sigma_tau | 1);
target += lkj_corr_cholesky_lpdf(L_phi | 4);
target += lkj_corr_cholesky_lpdf(L_tau | 4);
} // end model block
generated quantities{
// Corr. matrix from cholesky factors
// Varying individual effects
matrix[12, 12] c_ind;
c_ind = multiply_lower_tri_self_transpose(L_phi);
// Varying dyadic effects
matrix[12, 12] c_dyad;
c_dyad = multiply_lower_tri_self_transpose(L_tau);
// Avg. holding time per state (across dyads)
vector[3] avg_theta;
avg_theta[1] = mean(theta[, 1]);
avg_theta[2] = mean(theta[, 2]);
avg_theta[3] = mean(append_row(
to_vector(theta[, 3]), to_vector(theta[, 4])));
// Avg. transition probabilities (across dyads)
matrix[3, 3] avg_Gamma;
avg_Gamma[1, 2] = mean(Gamma[, 1, 2]);
avg_Gamma[1, 3] = mean(Gamma[, 1, 3]) + mean(Gamma[, 1, 4]);
avg_Gamma[2, 1] = mean(Gamma[, 2, 1]);
avg_Gamma[2, 3] = mean(Gamma[, 2, 3]) + mean(Gamma[, 2, 4]);
avg_Gamma[3, 1] = mean(append_row(
to_vector(Gamma[, 3, 1]), to_vector(Gamma[, 4, 1])));
avg_Gamma[3, 2] = mean(append_row(
to_vector(Gamma[, 3, 2]), to_vector(Gamma[, 4, 2])));
avg_Gamma[3, 3] = mean(append_row(
to_vector(Gamma[, 3, 4]), to_vector(Gamma[, 4, 3])));
}
To obtain prior samples, we create a dummy data set, that we fit to the prior model (these data do not matter, because the model has no likelihood function).
